{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Notes on Developing Fourier Series/Transforms for RCWA/PWEM\n",
    "PWEM and RCWA are dependent on expressing Maxwell's equations (in some number of dimensions) in k-space. Doing that either involves a fourier transform (aperiodic) or a fourier series (periodic system)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import signal\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy\n",
    "import cmath"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Heaviside Function Truncated to a Unit Cell\n",
    "The 1D fourier series for a step function is actually analytic, which we show here.\n",
    "It is extremely useful for many simulations of gratings in RCWA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEICAYAAAC3Y/QeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xm4ZFV57/Hvr4Yz9EQ3NspMq6BX\nMIraQhKNchUVjAIx4QbjVYwa9N7gkGgcYq4SjLma+2SORkkkMopcjUoiBjUGh6sMDUEFAWkR6IYG\nGnrmTDW894+1T1McztBn2DXs+n2e5zynqvauvdbetfZb715rVZUiAjMz6y+lTlfAzMzaz8HfzKwP\nOfibmfUhB38zsz7k4G9m1occ/M3M+lDPBH9JT5f0n5J2S3pHDtv/mqQzl3q7lkh6naSvd6DcF0i6\nQ9IeSae1u/w8SDpH0sWdrsekuc4dSZ+V9CdtrE9by+tVPRP8gfcCV0fEyoj4m6XeeEScHBEXLPV2\nu5WkkHRkTttel22/MvlYRFwSES/Po7w5nAv8XUSsiIgvd6D8wms9dyS9UdL3Ol2nbifpREk3SnpE\n0iZJ/61l2bGSbpA0kv0/tmWZJH1c0sPZ359J0kLq0EvB/wjglqXeaHYwF3wcWgPcYtbpNr1Y5xnk\n0m6WymLbXz9b6jbarjYv6WjgUuCDwH7AscAN2bIB4CvAxcAa4ALgK9njAGcBpwHPBp4FvAp464Iq\nEhFd/wd8C2gAY8Ae4GnZQbsQ2ArcDfwRUMrWPwe4uOX564AAKtn9q4GPAv8PGAWOzB57S8tz3gTc\nCmwHrgKOaFkWwO8CdwA/n6a+k+W9GbgH+E72+C8C3wd2AD8ETmh5zv7APwH3ZWV+uWXZ7wAbgW3A\nFcDBU+rytqwu24FPAMqWHQl8G9gJPAR8Pnv8O9nzHsmO528CJwCbgfcB9wMXAW8Evjdl3wI4Mrs9\nDPx5dvx3At/LHrsnW29P9vdLU7cF/DJwffa864Ffbll2NfCR7PXZDXwdWDtL+5j2+AA/A5rZa7wH\nGJzmue8D7s3KuR14acu+fTY7pj8B/gDYPN1xyO5/FviT7PYa4F9JbXN7dvvQKfs3tf3tB3wG2JLV\n50+A8gz7ew7wBeDzWb1vBJ7dsvwZWRk7SG98p2SPDwA3AW/P7pezOnxomjKenD1/8pz6R+DBluUX\nA+9q2Z+3ZOWOkc7VPcCOlmPzCeCrWX2vBZ46y+t5SlbvHdm2n9Gy7K7sNfsRMA5UgOdkx2B3dkwu\nm3wtsue8KtvvHaTz71lzbG/aNrGE8exS4CMzLHt5VrZaHrsHOCm7/X3grJZlbwauWVA9lnKn8vzj\n8cH5QtI75EpSsP0p8OaWk2Ou4H8PcEz2Yldbt096Z92YNeYK6Y3l+1NO/G+QAvbwNHWdLO9CYDkp\nkBwCPAy8knTF9bLs/gHZc76aNdw1WX1enD3+ElLgfi4wCPwt2ZtJS13+FVgNHE4KOJMN5XOk7KIE\nDAEvnPK81uB1AlAHPp6VM8zcwf8T2XE7hBRIfjl77mOOd7bu3m1lx2078Prs+L42u/+EltfnZ6Q3\n+eHs/sdmaBdzHZ+7gBNneO7TgU08+maxjiwoAR8DvpvV9TDgZvY9+D8B+HVgGal9/l8e+2Z+NY9v\nf18GPk1qL08ErgPeOkO9zwFqwG9kz30P8PPsdpXUdv+QFOxfQgpiT8+e+8zsWD8jaxvXMPObzD3A\n87LbtwN3kgXibNlzpp6bTN9mPkt6Yz4u299LgMtmKPNppKTkZdm+vDfbn4GW1/Om7DUZzvbxbuD3\nsvV/Izs2k6/Fc4EHgeNJbfTMbBuDM2xvxjYxTV3fT3pDmfZvllh2Jym5+THpzf5iYP9s2e8BX5uy\n/r8C785u7wSOb1m2Hti9oJi6VME5778pDaxMepc+umX5W0ljArBvwf/cWbb/NbI3kux+CRghy/6z\nbb1klrpOlveUlsfeB1w0Zb2rssZ4EClDXTPNtj4D/FnL/RVZ417XUpfWoH458P7s9oXAebRknS3r\nTRf8J4ChlsfeyAzBPzsmo7RknDMd76nbIgX966Y85wfAG1teiz9qWfY/gX+b4VjPdXzuYubgfyQp\nMJwIVKcsu5PsTTS7fxb7GPynKedYYPuUtnZuy/0nkdrzcMtjrwX+Y4btnUNLtpe9FluAX8n+7ifL\n2LPlnwPOabn/buA20pvAUbO044uA3wcOJAX/PyNdZU69KriauYP/P7bcfyVw2wxl/i/g8in7di/Z\nVXL2er6pZfmLSFfLrZny93k0+P89U7LsbF9ePMP2ZmwTS/VHOs/uIr3RrQC+CFzSsv+XTVn/ksnX\nj3RV9V9alh2VtUXNtx692te4lkff8SfdTcpA99WmWZYdAfy1pB2SdpCyFk3Z/mzPn26dI4DTJ7eZ\nbfeFpMB/GLAtIrZPs42DadnPiNhDumJorcv9LbdHSA0KUtYk4DpJt0h60xz13RoRY/uwX5BegyFS\nhj5fj9mnzNTXb6Z9mnVbMxyfaUXERuBdpGD6oKTLJB3cst3W129qfWckaZmkT0u6W9IuUjfbaknl\nltWmto0qsKWlbXyadAUwk73Pj4gmqcvu4Ml6Z4+11r31eFxAeoO+MiLumKWMb5OSghdl+3A18OLs\n77tTypjLQl/PJmlfZzr3DgbujSwSZlpfqyOAd0857w7Lnve47c3RJpbKKPBPEfHTrL3+KekNEVJ3\n2aop668iXb1Nt3wVsGfK/u+TXg3+D5GyuyNaHjuclCFAumxc1rLswGm2MdvB2kS65F7d8jccEd/f\nx+dPt84mUubfus3lEfGxbNn+klZPs437aNlPSctJ3Qr3TrPuYwuPuD8ificiDiZdGX1yjhk+U/fp\nMcdRUutxfIjUv/vUfdjOVI/Zp0zr6zcfCz4+ABFxaUS8MNtGkLq9IGXSh02pX6sRZm5j7yZ1Hxwf\nEatIwRPSG/HeoltubyJl/mtb2saqiDhmlqrvrVs2YHwo6VjcBxw2ZRB56rH9JKkr4RWSXjhLGd8m\nXUmckN3+HvACUvD/9gzPmXcQmmLq6ynSvrbWv7WMLcAhU2a8tL5Wm4CPTjnvlkXE52aq8yxt4jEk\n/WE2hXjav1n28UdTy2xxC/CsKfvzLB6dtHALabB30rNZ4ISGngz+EdEgdW98VNJKSUeQLk8n5z7f\nBLxI0uGS9gM+MM8iPgV8QNIxAJL2k3T6Iqt9MfBqSa+QVJY0JOkESYdGxBZSV9MnJa2RVJU0GTAu\nBX47m/41SMoSro2Iu+YqUNLpkg7N7m4nNbhGdv8B4ClzbOKHwDFZ2UOkbAjYm5GdD/yFpIOzffql\nrI5bSd1YM23/SuBpkn5LUkXSbwJHkwLSfC3m+Dxd0kuy542RMrLJ43M5qQ2syY7h26c8/Sbgt7L9\nPokUECetzLa1Q9L+wIdnq0f2+n8d+HNJqySVJD1V0otnedrzJL0mm6HyLtKbxzWkwdRHgPdm7egE\n4NWkQVAkvR54Hql75h3ABZKmzcKzq4JR4L+TxlF2kdrNrzNz8H8AOLRldsp8XQ78qqSXSqqS3kjH\nSV050/kBaazqHVlbeg1pbGHSPwBvk3R8NrNquaRflbRyuo3N0SYeIyL+NNIU4mn/ZtnHfyK12adI\nWkbqEp5s+1dn5b1D0qCks7PHv5X9vxD4fUmHZFck7yZ1q81bTwb/zNtJjfxOUkZyKSkYERHfIA2e\n/og0hWpeQSUivkR6t78su2y/GTh5MZWNiE3AqaSBuK2kjOQPePQ1eD3pauY2Up/ju7Ln/TupH/CL\npCznqcAZ+1js84FrsyzkCuCdEfHzbNk5pBN/h1rmGE+p809J8+S/SZpNNHX+9ntIg1bXk7rGPk7q\nBx4hm82Sbf8Xp2z3YdIMjHeTumjeC7wqIh7ax/1q3dZijs8gaWD3IVK3xBNJrw/AH5O6D35OCswX\nTXnuO0lBdQfwOtKA7aS/Ig0ePkQKyP+2D3V5A6kr8yekN+ovkLoEZ/IV0iytyYHz10RELSImSLNl\nTs7K/yTwhoi4TdLhWd3eEBF7IuJSYAPwl7OU823g4Yi4p+W+gP+cYf1vkTLR+yUt5PW8nfRm87dZ\n/V8NvDrbr+nWnwBeQ3oz2046Jv/csnwDaTbY32XLN2brzmS2NrEkIuJ8UhC/ltTGxklvxJP7cxqp\nPewgzTo8rWX/Pw38C+m8u5k0UeTTC6nH5JRAM5tFlkFfHBGHzrWuWS/o5czfzMwWyMHfzKwPudvH\nzKwPOfM3M+tDXfvlXWvXro1169Z1uhpmZj3lhhtueCgiDphrva4N/uvWrWPDhg2droaZWU+RtE+f\nRne3j5lZH3LwNzPrQw7+ZmZ9yMHfzKwPOfibmfUhB38zsz7k4G9m1occ/M3MushVt9zP1t3juZfj\n4G9m1iVGJuq89aIbOPP863Ivy8HfzKxLTNTTzyLfu2M097Ic/M3MukS9mb5luVLSHGsunoO/mVmX\nqDdS8C87+JuZ9Y9aI3X7VMv5h2YHfzOzLtFoOvM361q1RpP3fuGHbN4+0umqWMG4z9+si/3gZw9z\n+YbNfOCff9zpqljBTGb+lbKDv1nXKSmdmE3//rUtsck+/3LJff5mXWfyirzZ7Gw9rHh+99IbAaj2\nSuYv6XxJD0q6eYblkvQ3kjZK+pGk5y5FuWadoCzz/8GdD3PSX32HrbvHaTZ9FWALt3usxhvOv467\nH07jSO0Y8F2q3/D9LPB3wIUzLD8ZOCr7Ox74++y/Wc8JHg30t92/m+d/9JsAPPfw1Ry8epinrF3O\nssEKywbKPHHlIM1IU/cOWDlIMwIBT1g+yGitAcDywTLj2Sc7h6tlao0mQgwPlKk3m1TLJZoRRLad\nest0wFqjCYKBcomJehNJVMtiIttGtay9gaTeCIaqaZujEw2WD1ZoRBDNtE+lkojsakYl0m2lK51m\ndrtcEo1mIEFZot5yu5HtW0na2yU2uX7rbUlpm9lhLAkme9CU3ZZAaO+xnrwtRL3ZZPlAhR2jNQYq\nJZZVy+waq1EuiWakfvOBSjpO9Wba51q9SSPS7fFag2bA8ECZsVqDCBiqlva+HkPVMqMTDSQYLJcZ\nqdUpS1TKaZ1avUm1Utq7nWpZe59blnhkokFZQoIdIzVq2SXi5m0j7Blv0Izg1i272LJzjC07Rrlv\n59jj2lg7BnyXJPhHxHckrZtllVOBCyMigGskrZZ0UERsWYryzdqp1pg+y7/xnh3ceM+ONtfGiqiX\nMv+5HAJsarm/OXvsMcFf0lnAWQCHH354m6pmNj9vu+iGBT3vCcsHePiRCVYMVthvuMoDu8bYb7jK\n2hWDbNk5ytoVg6wcqrB19ziHrBmmXBK1RrBmWZVKqcRIrcGaZVUazWCi3mT/5QNM1JvUmmmd0YkG\njQhWDw8wWmsQEawaTo83I1g5VGWs1uCQ1cM0I9iSZZwDlRLVshiZaDBYKbNsoMx4vUGjma5EJhop\nwx2slPZ+90y1XGKi0aSkNDhZazQpZ91h9WYwUE5XBc1I2681mjSzLHysnrLt4Wp5bz2XDVQYmagD\nMDxQYWS8jpTWeWSigUgZ+chEgyeuGmR0osH2kQkO3G+IsYkGW/eMs9/wAGO1Bo1msHwwbW9yv3eN\n1hiZaLBm2QC7xmp7j9+2RyZoNIP9lw+wdfc4QbBm2QAP7BpDEquGKjywa5xKWQxXy/z8oUdYMZTC\n5r3bR1m9rEqtETy4a4zVywYYmaizc7TGYnsBK20Y8G1X8J/ubexxhycizgPOA1i/fr07Ua0rTV7i\nT3rri57C6esP45DVwwwPlDtUK+s1e8bTm9337niI796xlUuuvWfvMuWf+Lct+G8GDmu5fyhwX5vK\nNsvNreee5IBvC7JiMIXfk555ICc980COe/L+vPOym4BHv+MnT+2a6nkF8IZs1s8vAjvd32+97uVH\nP8mB35bMKc8+mEvekubBHHv46tzLW5LMX9LngBOAtZI2Ax8GqgAR8SngSuCVwEZgBPjtpSjXrJM+\nctozO10FKxBJvODItQyUS7Tj84NLNdvntXMsD+B3l6Iss25RakfHrPWfNjUrf8LXbIGWD7rLx/IR\nj58Ps+TaNeBrVhjPX7eGSqnEsgGfPrb0BNPMhVx6zvzNFsA9PpYXqS2x38HfbL78ZZ6WJyGiDY3M\nwd9sAZz5W17a1bYc/M3myYm/5a0dV5cO/mYLoHbNx7O+I9znb9aV2tEfa/1LkjN/s27lPn/LS8r8\nPeBrZtZfPOBr1p3c6WN5c7ePmVmfaVePooO/2Tx5vNfylAZ83edv1pXkEV/LiT/kZdalnPhb3jzP\n36xLOe+3vAgP+Jp1J3f6W44keZ6/Wbdyl7/lxZm/WZdy3m958oCvWRdz4m958oCvmVnf8Re7mXUl\nj/dantSmH/F18DdbAH/Iy/LiAV+zLtWOaXjWvzzga9bFnPdbnpz5m3Uh9/lbnoQ/5GXWtdzlb3mR\nnPmbdSVn/pYn/4C7WVdz6m/5aNdMMgd/s3ly4m95c7ePWZdyn7/lyQO+Zl2oHT+xZ/1Lber0d/A3\nWwAn/pYXyQO+ZmZ9R21KLZYk+Es6SdLtkjZKev80y98oaaukm7K/tyxFuWZmRdSOrsXKYjcgqQx8\nAngZsBm4XtIVEfGTKat+PiLOXmx5Zt3AA76Wl17q9jkO2BgRd0bEBHAZcOoSbNesK3m81/LUS9/q\neQiwqeX+5uyxqX5d0o8kfUHSYdNtSNJZkjZI2rB169YlqJpZPtrVL2v9J/2Ae/6WIvhPdxZMrfu/\nAOsi4lnAN4ELpttQRJwXEesjYv0BBxywBFUzW3r+SmfLU7vSiqUI/puB1kz+UOC+1hUi4uGIGM/u\n/gPwvCUo16xj3OdveWrHgO9SBP/rgaMkPVnSAHAGcEXrCpIOarl7CnDrEpRr1hHu87dctWnAd9Gz\nfSKiLuls4CqgDJwfEbdIOhfYEBFXAO+QdApQB7YBb1xsuWad5Mzf8tKmn/BdfPAHiIgrgSunPPah\nltsfAD6wFGWZdZoTf8uTv9XTrIt5to/lyV/sZmbWZ3ppnr9ZX/G3elqe/DOOZt3MvT6WE/+Au1mX\nct5veWrXTDIHf7MFcOJveXK3j1k3cupvOeuV7/Yx6zvtmott/UeSM3+zbuTE3/KU0goP+Jp1Jef9\nlhcP+Jp1Kc/zt7y528esS7nL3/LSSz/jaGZmS0SoZ77P36yvuNPH8uTM36yLudfH8tJLP+No1lc8\n3mt584CvWZfyh7wsN5K7fcy6UTu+cdH6V/o+fw/4mnUl5/2WF3/Iy6xLuc/f8uQBX7Nu5tTfcuQB\nX7Mu5Mzf8iT5l7zMupac+ltO/APuZmZ9yD/gbtbFPM3f8tKuq0oHfzOzLuM+f7Mu5O/zt1y528es\ne7nXx/Ii/K2eZl3Jeb/lyZ/wNetiHvC1XLnbx6z7uMvf8iT8IS+zruUPeVlePM/frEv5K50tT/4Z\nR7Mu5j5/y0tPfchL0kmSbpe0UdL7p1k+KOnz2fJrJa1binLNOsF9/pa3nvgxF0ll4BPAycDRwGsl\nHT1ltTcD2yPiSOAvgY8vtlyzTnLmb3nppW6f44CNEXFnREwAlwGnTlnnVOCC7PYXgJcqpx9B3TVW\n439cfANX3/5gHps3M8tdrwz4HgJsarm/OXts2nUiog7sBJ4wdUOSzpK0QdKGrVu3LqgyE/UmX7v5\nfu7ZNrKg55vNxb0+lqdD1yzjsP2X5V5OZQm2MV0GP/X82Jd1iIjzgPMA1q9fv6BzbLIg98tavtzv\nY/n436/5hbaUsxSZ/2bgsJb7hwL3zbSOpAqwH7BtCcp+nJx6k8z2cmJhRbAUwf964ChJT5Y0AJwB\nXDFlnSuAM7PbvwF8K3IezvY3L1qenGNYr1t0t09E1CWdDVwFlIHzI+IWSecCGyLiCuAzwEWSNpIy\n/jMWW+5M9nb75FWAmVuXFcBS9PkTEVcCV0557EMtt8eA05eirLlMZmRO/C1PTvyt1xXuE76Tn45z\n7Le8OLGwIihc8HdKZu3gPn/rdcUL/hkP+Fpe3LKsCAoX/J2RWTv4K52t1xUv+Gf/nfibmc2seMFf\nkwO+jv6WD3cpWhEUL/h3ugLWF9y9aL2ucMF/kpMzy4ublhVB4YL/3g95dbYaVnBO/K3XFS/4T37I\ny9HfcuK2ZUVQvOC/N/P3GWr58bfHWq8rXPA3y5tn+1gRFDb4+/w0M5tZ4YK/r8Ytb84rrAiKF/z3\nDvj6FLX8OMmwXle84O+T0vLmvMIKoHDBf5ITf8uTv9jNel3hgr9/xtHMbG7FC/7yh7wsX25aVgTF\nC/7Zf3/Iy/LksSXrdcUL/j4pLWeeSWZFULjgP8nnp+XJOYb1usIF/0d/zMUsH25bVgSFC/57OfW3\nHLl70XpdIYO/5OzM8uO8woqgmMG/0xWwwvNXOluvK2TwB2dnlh9PI7YiKGTwl+QT1HLlvN96XTGD\nP878zcxmU8zg7wFfy5ETCyuCYgZ/X5Rb3tzErMcVMviDszPLj5uWFUExg788I8Py5atL63WFDP4C\np2eWH7ctK4BFBX9J+0v6hqQ7sv9rZlivIemm7O+KxZS5b/Xy+Wn58me8rNctNvN/P/DvEXEU8O/Z\n/emMRsSx2d8piyxzTr4ktzy5S9GKYLHB/1Tgguz2BcBpi9zekvF3rluenF5Yr1ts8H9SRGwByP4/\ncYb1hiRtkHSNpBnfICSdla23YevWrQuulOTZPpYfty0rgspcK0j6JnDgNIs+OI9yDo+I+yQ9BfiW\npB9HxM+mrhQR5wHnAaxfv37Bp5hwn7/ly33+1uvmDP4RceJMyyQ9IOmgiNgi6SDgwRm2cV/2/05J\nVwPPAR4X/JeKJGdnZmazWGy3zxXAmdntM4GvTF1B0hpJg9nttcALgJ8sstxZOSmzPDmvsCJYbPD/\nGPAySXcAL8vuI2m9pH/M1nkGsEHSD4H/AD4WEbkGf/CMDMuXZ5RZr5uz22c2EfEw8NJpHt8AvCW7\n/X3gFxZTzrx5wNdy5JlkVgTF/YSvWY484Gu9rpjB32em5ch5vxVBIYM/+NLc8uX0wnpdIYO/v9vH\n8uS8woqgmMEfn6CWM3ctWo8rZvD3D7ibmc2qmMG/0xWwwnMbs15XyOAP7vaxfHgigRVFIYO/B3wt\nb+7yt15XyOAP/mI3M7PZFDL4yz/iazlxUmFFUczg3+kKWOH5i92s1xUy+IMzNMuHm5UVRSGDv3/G\n0fLmAV/rdcUM/vhDXpYPT/W0oihm8Hfmbzlz4m+9rpjBv9MVsMJyTmFFUcjgDz5JLV/u87deV8jg\nL/lDXpYPtysrikIGf/APuFu+/Gtx1usKGfwl3O9jZjaL4gZ/sxz4itKKopDBH5z4m5nNppDBX8gf\nxrFcuFlZURQz+Pv7/C1n7lq0XlfM4I8zNDOz2RQz+Dsts5z5K52t1xUy+IO7fSwfvqK0oihk8E/d\nPj5LLT++uLReV8jgjwd8LSee529FUcjg76TM8nL9XdsB2DFS63BNzBankMEfcOpvufjU1T8D4Ia7\nt3W4JmaLU8jgL/mXvCwfg9V0ykzUmx2uidniFDL41xpNNj64h52+NLcldvXtWwEYd/C3Hreo4C/p\ndEm3SGpKWj/LeidJul3SRknvX0yZ++Luh0f46QN7ePa5X8+7KOsj9cajAd+Zv/W6xWb+NwOvAb4z\n0wqSysAngJOBo4HXSjp6keXusz/68o+5/q5tbNo20q4irWC2PTLB5667hyM/+LW9jznzt15XWcyT\nI+JWmPMTtccBGyPizmzdy4BTgZ8spux9dfE193DxNfc87vHBSom1KwZpRrD/8gFWDlUYmWjwxJVD\n7L+8yu6xOkPVMquGKuwaqzNQLrHfsiq7x2qUJJ6wfIDtIzVWDlVYPlhh12iNwWqZZQNldozUWDZQ\nZrBSYudojRVDFQbKJXaP1Vk5lA756ESD/ZZVqTfS2MTwQJlKSawcqrJ9ZIKSYKIRrB6usmM0dV+t\nHKywZ7xOSWL5YJk943XKEsMDZR4Zb1ApicFqidGJBuWSqJZLjNcbVEolSiWo1YNqJf3KWb0ZDFZK\n1BtpdGSgUmKi3iQiGB4oM1ZLwW2oWmKs1kTAULXMaK3xmNslpeeO15pIMFgps3Kows7RdJxWDVe4\nf+cY4/Umw9UyI7UG1ZIolcToRINquYSUjke1UoIIRmsNBitlao3m3ueN1RuM15qszF6PRrPJfsNV\ntj2Sjs2q4Qrb9kww0UiPP7RnHIDlAxW27hmnJDFULfHgrnGq5RLLBys8uHsMgJLEz7buodkM9ozX\n2TVWn7Nd/ebzD1tMszTruEUF/310CLCp5f5m4PjpVpR0FnAWwOGHH55rpcbrTe7dMQrAlp1jLUt2\n5lquFcPbX3Jkp6tgtihzBn9J3wQOnGbRByPiK/tQxnSXBdNOxYmI84DzANavX5/LdJ0VgxUkeNKq\nIZ72pBVIYrBS4vD9lzFRb1Ipl1g5WGGi0aQksWKwzEQjUlZdKVFrNKmWS1TKJcpZp9lQpcxEo0ml\nlB6bzFabEYzXU7ZaawQT2e3xepNao8mqoSpjtQaNZrByqMrIRJ2DVw8jwX07xqhnZaVMvMFApYQQ\npZQgM1gpMdFoQkC1UqJWb4KgLFFvBuVSyvKbkerfaAbNaMnyCYaqadsRsGwgZfMR6UpkZKIBpMcf\nGa8jiaFKKWX/EgPlEmO1BqWSqJbFWK3BASuGGKs32DFS40mrBhmvN9n+yARrlg+we6yWlVNhZKKO\nBMPVdLtUEmWJkYk6jSasGKowMl6nEcHygXRV1ohgWXYV0GgGg5Uy9WaTiex415qR6qNUl2pZCLFn\nvM6ygTL1LLNfOVRhz1idnaM1KiWxa6zOg7vHWDmYrhju3jZCpSRuu3/3jO3I3x9lvW7O4B8RJy6y\njM1A6zXyocB9i9zmPvuVo9by12c8h/2XD7SryCVxzMH7dboKlrlvxyg/vncnb73ohk5XxWzJtGOq\n5/XAUZKeLGkAOAO4og3lctIxB3LRm4/vucBv3eXg1cO84pgDufXckzpdFbMls9ipnr8maTPwS8BX\nJV2VPX6wpCsBIqIOnA1cBdwKXB4Rtyyu2vvmU69/XjuKsT4xPFDudBXMlsxiZ/t8CfjSNI/fB7yy\n5f6VwJWLKcvMzJZOO2b7mBWaDpJvAAAD/0lEQVTG2f/1SG5/YOaBYLNe4eBvNg/vecXTO10FsyVR\nyO/2MTOz2Tn4m5n1IQd/M7M+5OBvZtaHHPzNzPpQIWf7/P3rnrv3F5fMzOzxChn8T/6FgzpdBTOz\nrub02MysDzn4m5n1IQd/M7M+5OBvZtaHHPzNzPqQg7+ZWR9y8Dcz60MO/mZmfUgR0ek6TEvSVuDu\nnDa/Fngop223g+vfeb2+D71ef+j9fcir/kdExAFzrdS1wT9PkjZExPpO12OhXP/O6/V96PX6Q+/v\nQ6fr724fM7M+5OBvZtaH+jX4n9fpCiyS6995vb4PvV5/6P196Gj9+7LP38ys3/Vr5m9m1tcc/M3M\n+lBfBn9JH5H0I0k3Sfq6pIM7Xaf5kvR/JN2W7ceXJK3udJ3mQ9Lpkm6R1JTUM9P1JJ0k6XZJGyW9\nv9P1mS9J50t6UNLNna7LQkg6TNJ/SLo1az/v7HSd5kvSkKTrJP0w24c/7kg9+rHPX9KqiNiV3X4H\ncHREvK3D1ZoXSS8HvhURdUkfB4iI93W4WvtM0jOAJvBp4D0RsaHDVZqTpDLwU+BlwGbgeuC1EfGT\njlZsHiS9CNgDXBgRz+x0feZL0kHAQRFxo6SVwA3AaT32GghYHhF7JFWB7wHvjIhr2lmPvsz8JwN/\nZjnQc++AEfH1iKhnd68BDu1kfeYrIm6NiNs7XY95Og7YGBF3RsQEcBlwaofrNC8R8R1gW6frsVAR\nsSUibsxu7wZuBQ7pbK3mJ5I92d1q9tf2GNSXwR9A0kclbQJeB3yo0/VZpDcBX+t0JfrAIcCmlvub\n6bHAUySS1gHPAa7tbE3mT1JZ0k3Ag8A3IqLt+1DY4C/pm5JunubvVICI+GBEHAZcApzd2dpOb659\nyNb5IFAn7UdX2Zf69xhN81jPXTUWgaQVwBeBd025ku8JEdGIiGNJV+zHSWp7F1yl3QW2S0ScuI+r\nXgp8FfhwjtVZkLn2QdKZwKuAl0YXDt7M4zXoFZuBw1ruHwrc16G69K2sn/yLwCUR8c+drs9iRMQO\nSVcDJwFtHYQvbOY/G0lHtdw9BbitU3VZKEknAe8DTomIkU7Xp09cDxwl6cmSBoAzgCs6XKe+kg2W\nfga4NSL+otP1WQhJB0zOzpM0DJxIB2JQv872+SLwdNJsk7uBt0XEvZ2t1fxI2ggMAg9nD13TSzOW\nJP0a8LfAAcAO4KaIeEVnazU3Sa8E/gooA+dHxEc7XKV5kfQ54ATS1wk/AHw4Ij7T0UrNg6QXAt8F\nfkw6fwH+MCKu7Fyt5kfSs4ALSG2oBFweEee2vR79GPzNzPpdX3b7mJn1Owd/M7M+5OBvZtaHHPzN\nzPqQg7+ZWR9y8Dcz60MO/mZmfej/A6xtEmCD6EPDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1f42e898358>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_order = 600;\n",
    "\n",
    "#determine the b_n coefs, which are analytic \n",
    "\n",
    "b_n_coefs = list();\n",
    "for i in range(1,n_order):\n",
    "    if(abs(i)%2 == 1):\n",
    "        b_n_coefs.append(1*4/(i*np.pi));\n",
    "    else:\n",
    "        b_n_coefs.append(0);\n",
    "\n",
    "#construct fourier series\n",
    "x = np.linspace(-np.pi, np.pi, 1000);\n",
    "L = np.pi;\n",
    "f_eval = list();\n",
    "for i in range(len(x)):\n",
    "    f_x = 0;\n",
    "    for j in range(1,n_order):\n",
    "        f_x +=b_n_coefs[j-1]*np.sin(j*np.pi*x[i]/L);\n",
    "    f_eval.append(f_x);\n",
    "\n",
    "plt.plot(x, f_eval);\n",
    "plt.title('fourier reconstruction of square box with orders = '+str(n_order))\n",
    "#plt.plot(x, np.sin(x));\n",
    "plt.show()\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fourier transform in python is really a fourier series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEICAYAAABS0fM3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAHF1JREFUeJzt3XucHGWd7/HPr7tnJneSkAlCEghI\nRAIri4wsLHpEQEVUYHdRQVZROYeX5+ARFZeLuguinpVzXMBFVhcMK+sFcJENEdllWS6u6AZJ0OV+\niVzMEEgmkBu5zEx3/84f9fRMzaQnM1PdSaf7+b5fDOl66qnqpy5d36qnqmfM3RERkXjlGt0AERFp\nLAWBiEjkFAQiIpFTEIiIRE5BICISOQWBiEjkFAQtxMzeZmZPNbodsTOziWb2UzPbYGb/1Oj2NAsz\nO8jMfmNmm8zs07WsRzM71sy6dzD+O2b2l6nh/2lmq83sNTPbs5blaEaFRjcgBmb2PLAXUEoVv8Hd\nV9Xzfdz9F8BB9ZxnrMzse0C3u38pw+SnkWzvPd29WNeGtbYLgPvc/XAAM/sIqfVoZh8D/ru7v7XW\nN3L3T1Zem1kbcAVwlLv/VyhzYIG7r6j1vZqBrgh2nfe7+5TUT11DwMxqCvVap5ch9gOeHikEtK5H\ntB/w2LDhEddjHe0FTBj23nFxd/3s5B/geeCEEcadTLIDrgfuAw5OjXPgwNTw94CvhtfHAt3AhcDL\nwPcrZan6+wA/AXqA54BPp8ZdCtwC/ADYSHKmNbxtE4G/AV4ANgD3AxPH0O7ngb8AHgY2A4tIPmz/\nAmwC/h2YEerOD8t5DrAKeAk4PzWvDuCqMG5VeN0xbB2cD6wJ03582LTfAH4PrAa+k2r/iNOGtvQD\nfcBrwE9D+YXAi2EZngKOr7LOvhym6w/Tng18DPglcCXwKvBVkpOwL4V1uwb4R2CPYevk48BKYB3w\nSeAtYZ2uB761g/3tUuDHYZ6bwnbqGud+dUFqvZwKnAQ8Hdr/hSr70c3hvR4CDgvj/gL4ybC2XQ1c\nVaXN95BcMW8L6+3GYevx3DCuFIbXh+lOAh4P7/0i8Pkx7hvfC9vhDST7qIf53gP8RxjeHMo+1Ohj\nyE4/RjW6ATH8MEIQpHbCdwJt4cO3AmgP40f7wBaBy0kOeBNJBQHJgWY58FdAO3AA8Czw7jD+0vAh\nOzXUnVilfdeQHOTnAHngj8N7jdbu54GlJAf/OeGD+BBweJj+HuCSUHd+WM4bgcnAH5AE1wlh/GVh\nXrOBTuBXwFeGrYPLQjtOArYwGDJXAUuAmcBU4KfAX49x2oF1HYYPIjko75Nq9+tH2N6XAj9IDX8s\nvNf/JumOnQh8IqyzA4ApwK3A94etk++QnKm+i+QguDish8o6ffsO3n9bWKY88NfA0tT4sexXfxXW\ny/8I2+NHYR0eEuZ9wLD96LRQ//MkJx1twN4k+8n0ULcQ2n3ECO2+j9QJyQjr8f5h07wEvC28ngG8\nebzbN7W+CyOto1b/UdfQrrPYzNaHn8Wh7EPAz9z9LnfvJzl7nUhywB2LMskBtdfdtw4b9xag090v\nc/c+d38WuA44PVXnP919sbuXh09vZjmSg9V57v6iu5fc/Vfu3jvGdl/t7qvd/UXgF8AD7v6bMP0/\nk4RC2pfdfbO7PwL8A3BGKD8TuMzd17h7D8kZ90dS0/WH8f3ufgfJGdxBZmYkB7HPuvur7r4J+D/D\nlr/qtCOs6xJJiC00szZ3f97dfzdC3WpWufvV7l4M6/pM4Ap3f9bdXwMuBk4f1m30FXff5u7/RnJA\nvTGsh8o6Hb4O0+539zvcvURytXjYONraD3wtbNubgFnAN919k7s/RnKF8aZU/eXufkuofwVJeB3l\n7i+RnF1/INQ7EVjr7svH0ZaxtHWhmU1z93Xu/tCwcWPdvlFTEOw6p7r79PBzaijbh6RrAAB3L5Oc\ndc4Z4zx73H3bCOP2A/ZJhc964AskZ+kVK3cw71kkH+hqB7uxtHt16vXWKsNThs0z3ZYXwnts917D\nxgG84kP7kLeEeXcCk4DlqeX/11A+2rTb8eSm4WdIzlLXmNlNZrZPtbojGL6uqy1XgaHbZ7zrMO3l\n1OstwIRx3Jt4JQRI5X2qtSX93gPLFvaFbga30Q3An4fXf04SSvX0ZyRn+y+Y2c/N7OjUuDFv39gp\nCBprFckBG4BwFjuPpK8Tkh13Uqr+64ZNv6NfHbsSeC4VPtPdfaq7nzTG6deSdAG8PkO7s5iXer1v\neI/t3mvYuB1ZS3LAOiS1/Hu4+1gPBNutG3f/kSdPrOwXxl8+xnlVm1+15Soy9IC7s4y2X43XwLYL\nV5JzGdxGi4E3mdmhwPuAH9bwPtW2yYPufgpJl9liknsjMk4Kgsb6MfBeMzs+PMJ2PtBL0g8O8Fvg\nw2aWN7MTgbePY96/Bjaa2YXheey8mR1qZm8Zy8ThzO564Aoz2ydMf7SZdYyh3Vn8pZlNMrNDSG6S\n3hzKbwS+ZGadZjaLpO/6B2Ns/3XAlWY2G8DM5pjZu8fYntUk/feEaQ8ys+PC8m8jCZnSSBOPwY3A\nZ81sfzObQtJtdbPvmsdNa9mvqjnCzP40XHF8hmRfWAoQrlhvIbnH8Gt3/30N77MamGtm7QBm1m5m\nZ5rZHqFbaiO1bZPh73XAqLVahIKggdz9KZLL5atJzmDfT/KYaV+ocl4oW0/Sp7y42nxGmHcpTPuH\nJDfv1gLfBfYYRxM/DzwCPEjytMjlQG4M7c7i5yQ3T+8GvhH6xSF5smMZydMyj5DcdP7qGOd5YZjn\nUjPbSPK00lj7iBeR9D1X7ul0AF8nWd6XSc5AvzDGeVVzPUk3yX+QbJ9tJDeTd4XM+9UIbiO5b7SO\n5P7Nn4YDc8UNJA8B1NotdA/J/YmXzWxtKPsI8HzYvp9ksBuqVpcCN4Tt/8E6zXO3Ze76wzTSOGY2\nn/CUyS46G5Y6MrNLSZ6uGfEAbGb7Ak8Cr3P3jbuqbTJ2uiIQkZ0m3DP4HHCTQmD3pW84ishOYWaT\nSfraXyB5dFR2U+oaEhGJnLqGREQi1xRdQ7NmzfL58+c3uhkiIk1l+fLla929c7R6TREE8+fPZ9my\nZY1uhohIUzGzF0avpa4hEZHoKQhERCKnIBARiZyCQEQkcgoCEZHIKQhERCKnIBARiZyCQKQGj3Rv\n4OHu9Y1uhkhNmuILZSK7q/d/634Anv/6exvcEpHsdEUgIhI5BYGISOQUBCIikVMQiIhETkEgIhI5\nBYGISOQUBCIikVMQiIhETkEgIhI5BYGISOQUBCIikVMQiIhETkEgIhI5BYGISOQUBCIikVMQiIhE\nTkEgIhI5BYGISOQUBCIikVMQiIhErm5BYGZ5M/uNmd0ehvc3swfM7Bkzu9nM2kN5RxheEcbPr1cb\nRERk/Op5RXAe8ERq+HLgSndfAKwDzg7lZwPr3P1A4MpQT0REGqQuQWBmc4H3At8NwwYcB9wSqtwA\nnBpenxKGCeOPD/VFRKQB6nVFcBVwAVAOw3sC6929GIa7gTnh9RxgJUAYvyHUH8LMzjGzZWa2rKen\np07NFBGR4WoOAjN7H7DG3Zeni6tU9TGMGyxwv9bdu9y9q7Ozs9ZmiojICAp1mMcxwMlmdhIwAZhG\ncoUw3cwK4ax/LrAq1O8G5gHdZlYA9gBerUM7REQkg5qvCNz9Ynef6+7zgdOBe9z9TOBe4LRQ7Szg\ntvB6SRgmjL/H3be7IhARkV1jZ36P4ELgc2a2guQewKJQvgjYM5R/DrhoJ7ZBRERGUY+uoQHufh9w\nX3j9LHBklTrbgA/U831FRCQ7fbNYRCRyCgIRkcgpCEREIqcgEBGJnIJARCRyCgIRkcgpCEREIqcg\nEBGJnIJARCRyCgIRkcgpCEREIqcgEBGJnIJARCRyCgIRkcgpCEREIqcgEBGJnIJARCRyCgIRkcgp\nCEREIqcgEBGJnIJARCRyCgIRkcgpCEREIqcgEBGJnIJARCRyCgIRkcgpCEREIqcgEBGJnIJARCRy\nCgIRkcgpCEREIqcgEBGJnIJARCRyCgIRkcgpCEREIqcgEBGJnIJARCRyCgIRkcgpCEREIqcgEBGJ\nXM1BYGbzzOxeM3vCzB4zs/NC+Uwzu8vMngn/zgjlZmZ/a2YrzOxhM3tzrW0QEZHs6nFFUATOd/eD\ngaOAc81sIXARcLe7LwDuDsMA7wEWhJ9zgG/XoQ0iIpJRzUHg7i+5+0Ph9SbgCWAOcApwQ6h2A3Bq\neH0K8I+eWApMN7O9a22HiIhkU9d7BGY2HzgceADYy91fgiQsgNmh2hxgZWqy7lA2fF7nmNkyM1vW\n09NTz2aKiEhK3YLAzKYAPwE+4+4bd1S1SplvV+B+rbt3uXtXZ2dnvZopIiLD1CUIzKyNJAR+6O63\nhuLVlS6f8O+aUN4NzEtNPhdYVY92iIjI+NXjqSEDFgFPuPsVqVFLgLPC67OA21LlHw1PDx0FbKh0\nIYmIyK5XqMM8jgE+AjxiZr8NZV8Avg782MzOBn4PfCCMuwM4CVgBbAE+Xoc2iIhIRjUHgbvfT/V+\nf4Djq9R34Nxa31dEROpD3ywWEYmcgkBEJHIKAhGRyCkIREQipyAQEYmcgkBEJHIKAhGRyCkIREQi\npyAQEYmcgkBEJHIKAhGRyCkIREQipyAQEYmcgkBEJHIKAhGRyCkIRDJas2lb1dcizUZBIJLRkV+7\nu+prkWajIBARiZyCQEQkcgoCEZHIKQhERCKnIBARiZyCQEQkcgoCEZHIKQhERCKnIBARiZyCQEQk\ncgoCEZHIKQhERCKnIBARiZyCQCSDG371/JjKRJqBgkAkg0uWPDamMpFmoCAQEYmcgkBEJHIKApFx\nKJedr9z++Ijjv3L745TKvgtbJFI7BYHIODy2aiOL7n9uxPGL7n+Ox1Zt2IUtEqmdgkBkjJ5bu5lT\nrrl/1HqnXvNLnu15bRe0SKQ+FAQiY7D02Vd4xzfuYyy9PmWH4/7m5/zn717Z+Q0TqYNCoxsgsrvq\n2dTLL1es5UuLH+W13uK4pz/juqVMbs/z1T85lLce2Enn1I6d0EqR2jUsCMzsROCbQB74rrt/vVFt\nkTi5O8Wy07Opl55Nvby6pY/nejaz/IV1/OyRl+ryHpv7Snz25v8aGH7vH+zNEfvNYP/Oycyc1E7n\n1A46p3ZQyBlmVpf3FBmvhgSBmeWBa4B3At3Ag2a2xN1HfhxDGs7dMTPcfVg5mCX/AniomzPDgVLZ\ncXygXt6MYjk5CLs7DrTlcjhOb3+Z/nKZchnyOaOQM/pKZbb2lXhlcx9t+eRgWSw7fcUyG7f2019y\n8jnoLZbpLZbZtK3Ilt4ivcUyfaVkeNO2fjZs7Wdzb5HVG3t5cf3WXbruKn72yEs7DJm5MyYye2oH\nkzsKTJvYxrQJbUydUKA9n6O9kGNKR4EpEwp0FHJ0FHKUytCWN6ZNbKO9kKOQS9ZPf8nZc3I7E9vz\ntOdzFMtO2R2zZF13tCW9wsWyY4BZsq4LOaPkg9vKMPI5wyBMn7yGwW0+PL8q+4iCrXk06orgSGCF\nuz8LYGY3AacAdQ2CDVv7+dSPHqrnLLdTKvvA44KVDw4kH5qyezKc/Ic7lFIH0crHJKk7dL6VgykM\nHmDL4QOaHu8OyaGUgbql8KF3oLe/jOOUyoNtSs+vv5QcdEvuhGMI/SU9/tgo3eu20r2uMSG1q1TC\nvOxJ2OcM2vI5cqngqJQnIWS0F3LkDHJmFPJJWWW/NwwL4yrSw5XidIhV5FL10nt9Jfw8jHCSE5tK\n3crnyyp1DfqKZToK+e2CsZp8ziiVnWLJKeSTsC07251kAby+cwqXnnzI6DOtQaOCYA6wMjXcDfxR\nuoKZnQOcA7DvvvtmexeHzRn6dseqv+Th4D/wdgOhkN4pkzNiQjnh4J08kw6DO6r74A42sAiVYR/c\nUSsHdGDgA1EuD55558OZ+Ja+EsVSOWlDCKWyJ8FVSgXCQDtrXyUio0qfaFQ+L9v6ywNXi+kDfd6S\nLrNcH0zuKGBAqc8HrlRyucGTLxgWAAP/S8oqVyk5GwwPGLyqqQSFGQMnZpUDfUXlqipXCZUwX3eY\n1J5nS9/ox5vKcaJyxdtbTD6TZjYkRCrt2tq38z+ZjQqCapk55LDk7tcC1wJ0dXVlOkXdY1Ibt/6v\nY7JMGp3KmUi1rp/BOoMbyd0Hzmp8YJwPfECSbp/BD04hl6NYLtNf8oEPTj6fnGH1Fcv0FcuU3AcO\nBMWSs62/RG8xBFnZ2VYsUSo7r20r0lsskc/lKJXL9JWcrX1FNm0r0lcq4w5b+ops6SuxYWs/vcUy\nG7b0s6m3yPotfaxav3W3veppyxt77zGRGZPbmTYh6R6aUMgzbWKBSe15JrUXMIP2fI6pEwpMbC/Q\nHtZjqexMaMszZUKBfM6YUMgPHMQ6CjkmtOUH1m1lXbcXki6nsjulyomNGW15G9hm7uHAZ1DIJScZ\nlavdyolQZV+w1EF4pDPj9D6m7qPdQ6OCoBuYlxqeC6xqUFuEoR/IkT6cQ4uTgUK+et1Cfvuy9pGe\nVm7QwzS9xRJbeku81lvk1c19vPDqFh5ftZGnV2/inifX7JT3PO6Ns3nDXlNZuM809ps5iZmT25nS\nUWBSR56OaiutwUbcZlWMtC9UowDYvTQqCB4EFpjZ/sCLwOnAhxvUFolURyE5+M6Y3M68mZM4bN50\nTj5sHyDptut5rZefP93DZT99PNPjowCT2/NccvIh/LcFncye2kEupwOg7H4aEgTuXjSzTwF3kjw+\ner2763f4ym4jlzP2mjaBD3bN44Nd87j7idWcfcOycc1j0VldHH/wXjuphSL1YyP1B+9Ourq6fNmy\n8X0IRept5atbeNv/vXdMdX9xwTuYN3PSTm6RyI6Z2XJ37xqtnn7FhMgYzZs5ies+Oupnius+2qUQ\nkKaiIBAZh3cu3GuHz4mbJXVEmomCQGScjn1D54jj3r6DcSK7KwWByDj93ZlHcMLBs7crf9uCWXz7\nzCMa0CKR2igIRMZpYnueAzqnbFf++s4pTGzf/b4LIDIaBYGISOQUBCJ1oi/LSrNSEIjUyfa/21Kk\nOSgIRDKo9kVMXRFIs1IQiGRwzIGztit7a5UykWagIBDJ4NiDZvOtDx8+MHz1GYfzjjdu/0ipSDNQ\nEIhkNLl98Hc2Tu7QY6PSvBQEIllZ+qVuEEjzUhCIZFTl7/SINCUFgUhGQ/6qWwPbIVIrBYFIRumD\nv/70ojQzBYFIRjbkHoFI81IQiGSUvkGsCwJpZgoCkYxMTw1Ji1AQiGQ09B5Bw5ohUjMFgUhWukcg\nLUJBIJKRKQmkRSgIRDLSPQJpFQoCkYx0j0BahYJAJCN9s1hahYJAJKMhXUO6JJAmpiAQyUhdQ9Iq\nFAQiGelXTEirUBCIZKZfMSGtQUEgkpHpDxJIi1AQiGSkewTSKhQEIhnp8VFpFQoCkYz0h2mkVSgI\nRDLSU0PSKhQEIhnpD9NIq1AQiGSkXzonrUJBICISOQWBSEZDf9dQ49ohUisFgYhI5GoKAjP7f2b2\npJk9bGb/bGbTU+MuNrMVZvaUmb07VX5iKFthZhfV8v4ijaT7AtIqar0iuAs41N3fBDwNXAxgZguB\n04FDgBOBvzOzvJnlgWuA9wALgTNCXZGmo64haRU1BYG7/5u7F8PgUmBueH0KcJO797r7c8AK4Mjw\ns8Ldn3X3PuCmUFdERBqknvcIPgH8S3g9B1iZGtcdykYq346ZnWNmy8xsWU9PTx2bKVIfenxUWkVh\ntApm9u/A66qM+qK73xbqfBEoAj+sTFalvlM9eLza+7r7tcC1AF1dXVXriDSSvlAmrWLUIHD3E3Y0\n3szOAt4HHO/ulQN2NzAvVW0usCq8HqlcREQaoNanhk4ELgROdvctqVFLgNPNrMPM9gcWAL8GHgQW\nmNn+ZtZOckN5SS1tEGkU3SyWVjHqFcEovgV0AHeF37641N0/6e6PmdmPgcdJuozOdfcSgJl9CrgT\nyAPXu/tjNbZBRERqUFMQuPuBOxj3NeBrVcrvAO6o5X1FdgdDfg21bhZLE9M3i0UyUteQtAoFgYhI\n5BQEIpnpT1VKa1AQiIhETkEgkpHuEUirUBCIZGQ7GBJpJgoCEZHIKQhEMjLT7xqS1qAgEBGJnIJA\nJCMb4bVIs1EQiGQ09KkhRYE0LwWBiEjkFAQiGZm+WSwtQkEgIhI5BYFIRvpmsbQKBYFIHejvEUgz\nUxCIiEROQSCSkbqGpFUoCEQy0ncHpFUoCEQyGvLNYmWCNDEFgUhG+maxtAoFgYhI5BQEIhnlU1cB\nOV0QSBMrNLoBIs2qc2oHH+qah+O8btqERjdHJDMFgUhGZsblp72p0c0QqZm6hkREIqcgEBGJnIJA\nRCRyCgIRkcgpCEREIqcgEBGJnIJARCRyCgIRkciZuze6DaMysx7ghRpmMQtYW6fmNAstc+uLbXlB\nyzxe+7l752iVmiIIamVmy9y9q9Ht2JW0zK0vtuUFLfPOoq4hEZHIKQhERCIXSxBc2+gGNICWufXF\ntrygZd4porhHICIiI4vlikBEREagIBARiVxLB4GZnWhmT5nZCjO7qNHtqRczm2dm95rZE2b2mJmd\nF8pnmtldZvZM+HdGKDcz+9uwHh42szc3dgmyM7O8mf3GzG4Pw/ub2QNhmW82s/ZQ3hGGV4Tx8xvZ\n7qzMbLqZ3WJmT4btfXSrb2cz+2zYrx81sxvNbEKrbWczu97M1pjZo6mycW9XMzsr1H/GzM7K2p6W\nDQIzywPXAO8BFgJnmNnCxraqborA+e5+MHAUcG5YtouAu919AXB3GIZkHSwIP+cA3971Ta6b84An\nUsOXA1eGZV4HnB3KzwbWufuBwJWhXjP6JvCv7v5G4DCSZW/Z7Wxmc4BPA13ufiiQB06n9bbz94AT\nh5WNa7ua2UzgEuCPgCOBSyrhMW7u3pI/wNHAnanhi4GLG92unbSstwHvBJ4C9g5lewNPhdd/D5yR\nqj9Qr5l+gLnhA3IccDtgJN+4LAzf5sCdwNHhdSHUs0YvwziXdxrw3PB2t/J2BuYAK4GZYbvdDry7\nFbczMB94NOt2Bc4A/j5VPqTeeH5a9oqAwR2qojuUtZRwKXw48ACwl7u/BBD+nR2qtcq6uAq4ACiH\n4T2B9e5eDMPp5RpY5jB+Q6jfTA4AeoB/CN1h3zWzybTwdnb3F4FvAL8HXiLZbstp7e1cMd7tWrft\n3cpBYFXKWupZWTObAvwE+Iy7b9xR1SplTbUuzOx9wBp3X54urlLVxzCuWRSANwPfdvfDgc0MdhdU\n0/TLHLo2TgH2B/YBJpN0jQzXStt5NCMtY92WvZWDoBuYlxqeC6xqUFvqzszaSELgh+5+ayhebWZ7\nh/F7A2tCeSusi2OAk83seeAmku6hq4DpZlYIddLLNbDMYfwewKu7ssF10A10u/sDYfgWkmBo5e18\nAvCcu/e4ez9wK/DHtPZ2rhjvdq3b9m7lIHgQWBCeNmgnueG0pMFtqgszM2AR8IS7X5EatQSoPDlw\nFsm9g0r5R8PTB0cBGyqXoM3C3S9297nuPp9kW97j7mcC9wKnhWrDl7myLk4L9ZvqTNHdXwZWmtlB\noeh44HFaeDuTdAkdZWaTwn5eWeaW3c4p492udwLvMrMZ4UrqXaFs/Bp9w2Qn34w5CXga+B3wxUa3\np47L9VaSS8CHgd+Gn5NI+kbvBp4J/84M9Y3kCarfAY+QPJHR8OWoYfmPBW4Prw8Afg2sAP4J6Ajl\nE8LwijD+gEa3O+Oy/iGwLGzrxcCMVt/OwJeBJ4FHge8DHa22nYEbSe6B9JOc2Z+dZbsCnwjLvgL4\neNb26FdMiIhErpW7hkREZAwUBCIikVMQiIhETkEgIhI5BYGISOQUBCIikVMQiIhE7v8De6+WXAG1\nx/4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1f42fcc52e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xt8VNW5//HPMzO5kRu3gEgCwTaI\nYLgEBEREeCGKN2hRqxyrYGupWD3+2opaeyqKv/68tt6O1FIV0GNVlKq0peLxHC3QqhAQsQIB5BoQ\nCQmXhNwmM8/vj4SYTAIZ4gx7svO8Xy9eZM9amflOMvNkzd5rry2qijHGGHfxOB3AGGNM5FlxN8YY\nF7LibowxLmTF3RhjXMiKuzHGuJAVd2OMcSEr7sYY40JW3I0xxoWsuBtjjAv5nHrgrl27anZ2tlMP\nb4wxbdKaNWsOqGpGS/0cK+7Z2dnk5+c79fDGGNMmicjOcPrZbhljjHEhK+7GGONCVtyNMcaFrLgb\nY4wLWXE3xhgXarG4i8gLIrJfRP51nHYRkadEZKuIrBeRvMjHNMYYczLCGbkvACaeoP0SIKfu3wzg\nd988ljHGmG+ixXnuqrpcRLJP0GUy8KLWXq/vIxHpKCI9VPXLCGU030DF0VLWLX4EqT5KaUJ3Puv+\nHQAGf7mIDv6DjfoeSuzJhm6XA5C392USa8oatRd36ENB14sAOKdwAXHBqkbt+5P7srXLOABG7v4D\nHg02av8yZQDbO49GtIZzdz/fJGth2mB2dRyBL1DJ8D0Lm7TvTD+HPel5JNQcYejeV5q0b+80ii9T\nc+lQXczgfW80ad/SeSxFKWeSWrWP3K/ebtJe0PVCijt8i44Vu+lftLRJ+4aMSzmUlEWX8i8488B7\nTdo/6z6Z0oTTyCgrIKfkgybt6067ivL4LvQo/Yw+B//ZpH3N6VOp8qXR8/Baeh9e3aR9Vc9p1HgT\n6XXoYzKPrGvS/mHWD1Hx0adkJT3KPm/UFhQPH2X9CIBvF79Pt6ObG7X7PQmszpwOwJkH3qVL+fZG\n7ZW+FNaefh0A/ff/hY6Ve75uTEgh73v3EJ+Q2CSTcU4kTmLqCexusF1Yd1uT4i4iM6gd3dOrV68I\nPLRpyZZVyzh321MA5Af78vSGswH4W9xr9JXCRn1X6tk8/Xk/AJbH/5GeHGjUvix4Dk/XnAHAmvgX\n6cjRRu2Lg+fzdE0WAD+Jn08cgUbtLwYu4unAacRRw+3xLzTJuiYwiacDXehIKbc2075iRylzA2n0\npIiZzbS/s72GhcEkzpRdzIhr2v7Gdh+Lg17yZDM/bKZ94fZU3gkq58t6pjfTPndbBiu0iomeVVzv\na9r+6LYs1mpfrvSs5Lpm2u/blkOB9mKaZznX+F5s0n7HFwPZQwa3eP+Xq72vN2mf+cVwDpHKHd73\nmOJd0qT9hi8uwI+P2d5lTPa+26jNj5frttb+4X3E91eGe1Y0aj9EMjdsGV37PH1LGO5p/MdlD125\nsWAEAAvj3mR43V5aj9Reg3nTp+fTb/iEJpmMcyScC2TXjdz/oqpnN9P2V+BBVV1Zt/0/wJ2quuZE\n9zls2DC1M1Sj79P/fZVBy3/M5klL6Jt3gdNxjMusXbOKz958hMFX38OgQXa47VQQkTWqOqylfpGY\nLVMIZDXYzgT2RuB+TQSE8bfbmFarTD+D2TU3Upna2+koJkQkivsS4Ia6WTMjgcO2vz12lKeewSP+\n7+HvcJrTUYwbaZA4aiDk+IpxXjhTIV8BPgTOFJFCEfmhiNwsIjfXdVkKbAO2An8AbolaWnPSKlJ7\nMzfwHQIp3Z2OYlwobf/HbEm8gdSvmh4ANs4KZ7bM1BbaFfhJxBKZiBJ/BT0pQgJVLXc2ppVs71/s\nsTNUXa7z/g/5R+LtJB7c3HJnY1pJrLzHHCvubmdHVE0USV0JsZdZ7LHi3m6I0wGMMaeQFXfXqx1S\nidV2EwXVKaczt2YSVck9nI5iQlhxd7mvT1Kz6m4iryq1F4/UXEtlqp1xHmusuLtcWfqZ3Oe/AX+y\nzXM3kSdBPx0pRYJ+p6OYEFbcXa48JYsFgYkEkro6HcW4UGrRJ6xL/LHNc49BVtxdzusvo6/sRmoq\nnY5ijDmFrLi7XJf9H/Juwl0kHNnmdBTjQseO5IjNhYw5Vtxdzt5zJpq0bhqWvcxijxV3l7MzB41p\nn6y4u17dPHf7VZso8Kdm8oj/e1Sk2FTIWGPveGNMq1Unn87cwHeoSslqubM5pay4u9yRjv2Z5Z9B\nTYrNczeR56mpoidFeGw2Vsyx4u5y5R0yeT0wlmBiJ6ejGBdKLv6MfyTeTmqRzXOPNVbcXS6u6iB5\nshnxVzgdxRhzCllxd7nORR/zp4T7iC/b6XQU40bHJrrbnNuYY8XdGPMN2IJ0scqKezsh9iY0pl2x\n4u5yX39atuJuIq8qtRez/dM4mnaG01FMCCvuxphWq+nQjYWBi6lK7ul0FBPCirvLHeo8mJnVtxNI\nsSvlmMjzBCroK7vx+sucjmJCWHF3uYoO3flbcATBhDSnoxgXSirZyLsJd5G2f43TUUwIK+4ul1BZ\nxPme9Xj8R52OYow5hay4u1znotW8FP8QvrK9TkcxLlQ/zd1WH405VtyNMd+AzcKKVVbcXe/Ykr/G\nRJ6IvbJilRV3l/t6mru9CU3kVaX2YpZ/BkfT+zodxYSw4t5uWHE3kRdI6sLrgbFUdbCptrHGirvL\nHehyDjdU30WNzXM3UeD1HyVPNuOtPuR0FBMirOIuIhNFpEBEtorI3c209xKR90XkExFZLyKXRj6q\naY2qxAyWBwehcclORzEulHBoK39KuI/0A2udjmJCtFjcRcQLPANcAvQHpopI/5Bu/wEsUtUhwLXA\n3EgHNa2TWP4lF3tW4bEzCE0UyNdzIU2MCWfkPhzYqqrbVLUaeBWYHNJHgWOnQKYDNqk6RnQpXsPv\n45/AV/6V01GMK9VWd6vtsccXRp+ewO4G24XAiJA+9wHvishtQDJwYUTSmYixw6kmmsQu1hFzwhm5\nN1cXQn+TU4EFqpoJXAq8JCJN7ltEZohIvojkFxUVnXxac/LsTWeiSWzkHqvCKe6FQFaD7Uya7nb5\nIbAIQFU/BBKBrqF3pKrzVHWYqg7LyMhoXWLTKmpjdxMF1am9mFl9O0c6DXA6igkRTnFfDeSISB8R\niaf2gOmSkD67gPEAInIWtcXdhuYxxM4kNNGgiR35W3AE1UndnY5iQrRY3FW1BrgVWAZspHZWzOci\nMkdEJtV1+znwIxH5FHgFmK5q+wNiwVfdRjGl6j5bz91Ehae6lPM96/FVFjsdxYQI54AqqroUWBpy\n270Nvt4AnBfZaCYSqhM6s1b7gi/J6SjGheKP7OSl+IdYV9yL2pnSJlbYGaou16FsF1d6liPVpU5H\nMS5mH9RjjxV3l+tcso7fxD+Lt8I+NhvTnlhxbyfseKqJCnthxSwr7m5nH5dNFFlpj11W3NsN+1Wb\nyKtO68UN1XdxqMtgp6OYEPaOdzkbt5to0vg0lgcHUZXQ5JxF4zAr7i73ZfexTKx6iEDKaU5HMS7k\nqT7CRM8qEipsYbpYY8Xd5fzx6WzSXogvwekoxoXiSvfwbPwTdCz51OkoJoQVd5dLObKVG7zLkCqb\n524i7+vJMrYDMNZYcXe5zoc+Y07cQjxVdhk0E0U2KyvmWHE3xrTasZG7lfbYY8Xd9extZ0x7ZMW9\nnbATCU00VKf1ZkrVfZR0He50FBPCirvb2cDdRFNcMmu1L9UJnZxOYkJYcXe5wh4XMbrqCVvP3USF\np+owV3qWk3S00OkoJoQVd5fz+zpQqN0Qb5zTUYwLxZXv4zfxz5JW8i+no5gQVtxdruORAm7xvm3z\n3I1pZ6y4u1ynQ59zZ9xreKqPOB3FuJjYwZ2YY8Xd9exNZ6JH6hb9VXudxRwr7u2GzYU0pj2x4u52\ndaeF2zx3Ew016b24uOoh9mec53QUE8KKe3th1d1EgfoSKdBe1MSnOR3FhLDi7nLbe17BkMpn0Q7d\nnY5iXMhTeYgbvMtILt3udBQTwoq7ywU8CRwkDTxep6MYF/JWFDMnbiFphzc4HcWEsOLucl0OfcYs\n36s2z92YdsaKu8t1PLKJn/iW4Kk56nQU42a2nnvMseLuevamM9Fjx+ljlxX3dkJsnrsx7YoVd7ez\nj8smimrSejO66gn2nTbW6SgmhBX39sI+P5to8MZRqN0I+JKdTmJCWHF3uS1ZV9K3ciHB5AynoxgX\n8lQe5Bbv26Qc3uJ0FBMirOIuIhNFpEBEtorI3cfp8z0R2SAin4vIHyMb07SWipdq4mzkbqLCW3WI\nO+NeI/3wJqejmBC+ljqIiBd4BpgAFAKrRWSJqm5o0CcH+AVwnqoeFJFu0QpsTk63kjXM9r2Kp3ok\nkOh0HGPMKRLOyH04sFVVt6lqNfAqMDmkz4+AZ1T1IICq7o9sTNNaHUu3cKNvGdRUOB3FuNDXnwft\nwH2sCae49wR2N9gurLutob5AXxH5h4h8JCITm7sjEZkhIvkikl9UVNS6xMaYGGK7+2JVOMW9ud9e\n6J9pH5ADjAWmAs+JSMcm36Q6T1WHqeqwjAw7wHdq1C3563AKY8ypFU5xLwSyGmxnAnub6fO2qvpV\ndTtQQG2xN7HCDqiaKAik92JI5bMU9rjI6SgmRDjFfTWQIyJ9RCQeuBZYEtLnLWAcgIh0pXY3zbZI\nBjWto3ioVq+N3E10eLwcJI2gN8HpJCZEi8VdVWuAW4FlwEZgkap+LiJzRGRSXbdlQLGIbADeB2ap\nanG0QpvwFWR9j75VL6EdujodxbiQt/IQs3yvkn7QlvyNNS1OhQRQ1aXA0pDb7m3wtQI/q/tnjGkn\npOoIP/EtYXXpSMB2zcQSO0PV5XoUf8SjvmcRvy35a0x7YsXd5dLKtnG1bzkSqHY6inGhr4/T2zz3\nWGPFvZ2wA6omOuyVFausuLuejahM9NnK0rHHinu7YSMsE3mBtCz6Vi5kZ+blTkcxIay4u1zQE88h\nTUY8VtxNFIjUrTrqdTqJCWHF3eU29byKwVV/QBObrAZhzDfmqTrIfb4FdDr4qdNRTAgr7u2EXUPV\nRIOnuozpvndJK7MT0mONFXeXyzywnLlxT0B1mdNRjIuJHVGNOWGdoWrarrTyXYz0ruJwMOB0FONC\nIjY+jFX2m2kvbK+MiSYbucccK+5uZ286E01Su+qo2pLSMceKezsh9uYzURBMPZ2+VS+xPeu7Tkcx\nIay4u5zf24FC7WoX6zBRZR8QY48Vd5cr6DmF0VVPQUKq01GMC3mqDvKo71kyilc7HcWEsOLeTti4\n3USD+Cu42rec1PJdTkcxIay4u1zvr97jxbgHwV/udBTjQvUnx9l+mZhjxd3lUiv3Msb7GaJBp6MY\nN7JjOTHLirvL2XjKnApqI/eYY8Xd5cTKu4ki8Xg5qCkEPfFORzEhrLi3EzbP3URDMLkbQ6rmsT1z\nstNRTAgr7i5X6UunIJgJtgaIiYJjQwa1T4gxx97xLlfQYzIXVz8CcR2cjmJcSCoPMTfuCU4r+qfT\nUUwIK+4ud2w8ZXtlTDRIoJpLvatIKS90OooJYcXd5XL2LWVx/GzwVzgdxbiRDRpilhV3l0uu2s9Q\nzxZ7D5oos33uscaKu8vZW85Elw0bYpUVd2NM64mXQu2K32sH7GONFXfXqx272zx3Ew3aoQujq55i\n++mXOx3FhLDi7nJH4zNYG/y2TZcxUWEvq9gVVnEXkYkiUiAiW0Xk7hP0u0pEVESGRS6i+SY2n3YZ\nU6rngC/R6SjGhaTyMC/GPUjm/g+cjmJCtFjcRcQLPANcAvQHpopI/2b6pQL/Dnwc6ZCm9Ww9JxNV\nQT9jvJ/RoXKf00lMiHBG7sOBraq6TVWrgVeB5haSeAB4BKiMYD7zDZ315Zu8E38X1FQ5HcW4kO2V\niV3hFPeewO4G24V1t9UTkSFAlqr+JYLZTAQk+Q/Sz7Pb9o2a6LKPiDEnnOLeXFmo/02KiAd4HPh5\ni3ckMkNE8kUkv6ioKPyUptXE3nMmisQWpItZ4fxmCoGsBtuZwN4G26nA2cAHIrIDGAksae6gqqrO\nU9VhqjosIyOj9anNSTg2FdLehCYKPF4KgplUxaU7ncSECOcdvxrIEZE+IhIPXAssOdaoqodVtauq\nZqtqNvARMElV86OS2JyUr5ditf0yJvI0sSMXVz/CttMucTqKCdFicVfVGuBWYBmwEVikqp+LyBwR\nmRTtgOabKU3syfJArk1INlFhr6rY5Qunk6ouBZaG3HbvcfqO/eaxTKQUZEzkN/4z2OIN61dtzMmp\nKmVx/GyO7LsJuNnpNKYB2xFrjGk1j9Yw1LOFpOpip6OYEFbcXW7g3tdYEX87EqxxOoox5hSy4u5y\n8TWlZHls2qmJMpvnHnOsuLcbdujLRIFNsY1Z9ptxOVFb8tdEkXhZG/w2FfFdnU5iQlhxd7kGpxI7\nGcO4VWIqU6rnsO20i51OYkJYcXe5g0nZ/C1wju2UMaadseLucpszLmSm/6e2b9RER1Up78Tfxbf2\nLmm5rzml7B1vjGk1Qenn2U2i/7DTUUwIK+4uN7TwRdYmzECwqWommuz1FWusuLucL1BJZymzA6om\nKmwWVuyy4u5ygk2FNKY9stWkXM4+LJuo8nhZHsilOuE0p5OYEDZyN8a0msQnc4P/F2zvNsHpKCaE\nFXeXK0rOYXHgfKdjGGNOMSvuLre18zh+7p/pdAzjVtVHWRF/O2fuecPpJCaEFfd2wI6lmmgRlCxP\nEfGBo05HMSGsuLvcyN1/YGP8NKdjGLezJX9jjhV3lxMNEIddqMNEh02xjV1W3F3Ozkw1pn2yee4u\npwpqa0KaaPF4a1cdTerldBITwkbuxphWE18SM/0/ZUfGOKejmBBW3F3uy5QBvBy0E0xMdNnx1Nhj\nxd3lvuh8PnMC052OYVxKaipYmzCD3MJXnI5iQlhxdznRGuJttoyJos5ShjdY5XQME8KKu8udt2se\nn8ZNdzqGcatjUyE16GwO04QVd2OMcSEr7q5nR7pMNNk021hl89zbBXsDmugQj4fFgfNJSP6W01FM\nCBu5G2NazxvHz/0z2dlljNNJTAgbubvcrtQ8PgmWcovTQYwxp1RYI3cRmSgiBSKyVUTubqb9ZyKy\nQUTWi8j/iEjvyEc1rbG947k8odc4HcO4lAT8FCTcwJBdC5yOYkK0WNxFxAs8A1wC9Aemikj/kG6f\nAMNUdSDwBvBIpIOa1vEFykmnzOkYxsUSpAZRO5ci1oQzch8ObFXVbapaDbwKTG7YQVXfV9Xyus2P\ngMzIxjStdd6e51nhtSsxmeiwJX9jVzjFvSewu8F2Yd1tx/ND4G/NNYjIDBHJF5H8oqKi8FMaY2Ka\n2OIyMSec4t7cn+Zmf5Mi8n1gGPBoc+2qOk9Vh6nqsIyMjPBTmlazN52JqrqRu73KYk84s2UKgawG\n25nA3tBOInIh8EvgAlW1hSZiiK3nbqJFEBbWTKBjSuhhOOO0cEbuq4EcEekjIvHAtcCShh1EZAjw\ne2CSqu6PfEzTWjaiMlHl8TK75kZ2dh7ldBITosXirqo1wK3AMmAjsEhVPxeROSIyqa7bo0AK8LqI\nrBORJce5O3OKbet4Lv+pVzkdw7iVKnHUIBpwOokJEdZJTKq6FFgactu9Db6+MMK5TITsTBvGS5rB\nLKeDGFcSgS2JN/DRrh9jM6Bjiy0/4HIJ/sP04IDTMYzL2e6/2GPF3eXO27uAP3t+5nQM41JS/7+V\n91hjxd3lFFsT0kSTTYWMVVbcXc4KuzHtk60K2Q7YPHcTLSLwu5or6JY20OkoJoSN3F3OPi6bqBLh\n4Zqp7Eof4XQSE8JG7i5X0HEM736ZzOy6bb/fT2FhIZWVlY7mMm1bYmIimZmZxMXF0ZFSfAF7PcUa\nK+4utytlMK9Ll/riXlhYSGpqKtnZ2bain2kVVaW4uJjCwkL69OnDusQf8+Gem4DfOB3NNGC7ZVwu\nxX+AMyis366srKRLly5W2E2riQhdunSxT38xzoq7y43a91+8zC8b3WaF3XxToa8hW3009lhxNzHr\n0ksv5dChQyfsc++99/Lee++16v4/+OADLr/88lZ9r/laUG2wEItsn7vrtb0RlaqiqixdurTFvnPm\nzDkFiYxpe2zk7nYKGoO7YX77299y9tlnc/bZZ/PEE0+wY8cOzjrrLG655Rby8vLYvXs32dnZHDhQ\nuy7OAw88QL9+/ZgwYQJTp07lscceA2D69Om88cYbAGRnZzN79mzy8vLIzc1l06ZNAKxatYpRo0Yx\nZMgQRo0aRUFBgTNP2qV+E7ianennOB3DhLCRezt2/58/Z8PeIxG9z/6npzH7igEn7LNmzRrmz5/P\nxx9/jKoyYsQILrjgAgoKCpg/fz5z585t1D8/P5/FixfzySefUFNTQ15eHkOHDm32vrt27cratWuZ\nO3cujz32GM899xz9+vVj+fLl+Hw+3nvvPe655x4WL14csefc3s0NfIdb07/tdAwTwoq7y33WeQJv\nFXXnIaeDNLBy5Uq++93vkpycDMCUKVNYsWIFvXv3ZuTIkc32nzx5MklJSQBcccUVx73vKVOmADB0\n6FD+9Kc/AXD48GGmTZvGli1bEBH8fn+kn1K71pMiEmq6Ox3DhLDi7nJ7kgewlPRmi3tLI+xo0ePM\nrDhW7MPt35yEhAQAvF4vNTU1APzqV79i3LhxvPnmm+zYsYOxY8eeXGBzQh/E/5TVe6ZRe/lkEyts\nn7vLdazaywC2Oh2jkTFjxvDWW29RXl7O0aNHefPNNzn//POP23/06NH8+c9/prKykrKyMv7617+e\n1OMdPnyYnj17ArBgwYJvEt2YNsOKu8udu/9VnuXXTsdoJC8vj+nTpzN8+HBGjBjBTTfdRKdOnY7b\n/5xzzmHSpEkMGjSIKVOmMGzYMNLT08N+vDvvvJNf/OIXnHfeeQQCdjm4aLD13GOPnMxH3kgaNmyY\n5ufnO/LY7cnH//kD+h14h/T79gKwceNGzjrrLIdTnbyysjJSUlIoLy9nzJgxzJs3j7y8PKdjtWvH\nXkv+2Z3Jz7yec3/0pNOR2gURWaOqLe4Ds33u7YAblvydMWMGGzZsoLKykmnTpllhN6YFVtxdzx0f\nl//4xz86HcEcx/8LXM+3Oo7mXKeDmEZsn3s74IaRu4ldLwUvZm9qrtMxTAgr7i63pssV3C8znY5h\nXCxHdpNcXex0DBPCirvLfZnUl7+LnRpuoudt3z3k7XvV6RgmhBV3l8uo3MEwPnc6hjHmFLPi7nIj\nil7n4eBvnY4RcxYsWMDevXsjdn9PPPEE5eXlJ/U9J1pyeOrUqQwcOJDHH388EvGOa8WKFQwYMIDB\ngwdTUVHBrFmzGDBgALNmzeKtt95iw4YN4d2RO47bu4oVd+MYVSUYDDry2Ccq7q050ak1xf149u3b\nxz//+U/Wr1/PT3/600Ztx5ZUiJSXX36ZO+64g3Xr1pGUlMTvf/971q5dy6OPPhp2cbcD9rHJirs5\npZpb2vfdd9/l3HPPJS8vj6uvvpqysjIAVq9ezahRoxg0aBDDhw+ntLSUyspKbrzxRnJzcxkyZAjv\nv/8+UFusp0yZwsSJE8nJyeHOO+8Eagv19OnTOfvss8nNzeXxxx/njTfeID8/n+uuu65+xJqdnc2c\nOXMYPXo0r7/+OmPHjuXYSXYHDhwgOzu7/v7uuOMOcnNzGThwIE8//TRPPfUUe/fuZdy4cYwbNw7g\nuM/pnXfeoV+/fowePbp+YbNQF110Efv372fw4MGsWLGCsWPHcs8993DBBRfw5JNPsnPnTsaPH8/A\ngQMZP348u3btAmqXP545cybjxo3jjDPO4O9//zs/+MEPOOuss5g+fXqTx3nuuedYtGgRc+bM4brr\nrmPSpEkcPXqUESNGcP/997NkyRJmzZrF4MGD+eKLLyLzAjCnzrELI5zqf0OHDlUTfR89dYMemJ1V\nv71hw4bGHV64tOm/j+fVtlUdbb597X/VtpcdaNrWgu3bt6uI6IcffqiqqkVFRXr++edrWVmZqqo+\n9NBDev/992tVVZX26dNHV61apaqqhw8fVr/fr4899phOnz5dVVU3btyoWVlZWlFRofPnz9c+ffro\noUOHtKKiQnv16qW7du3S/Px8vfDCC+sf/+DBg6qqesEFF+jq1avrb+/du7c+/PDD9dsN24uKirR3\n796qqjp37lydMmWK+v1+VVUtLi6u//6ioqITPqeKigrNzMzUzZs3azAY1Kuvvlovu+yyZn9GAwYM\naJRl5syZ9duXX365LliwQFVVn3/+eZ08ebKqqk6bNk2vueYaDQaD+tZbb2lqaqquX79eA4GA5uXl\n6SeffNLksaZNm6avv/56/XZycvJx20Idey3N+o+79YVFi4/bz0QWkK9h1FgbuZtTruHSvh999BEb\nNmzgvPPOY/DgwSxcuJCdO3dSUFBAjx49OOec2pk+aWlp+Hw+Vq5cyfXXXw9Av3796N27N5s3bwZg\n/PjxpKenk5iYSP/+/dm5cydnnHEG27Zt47bbbuOdd94hLS3tuLmuueaaFrO/99573Hzzzfh8tef/\nde7cuUmf4z2nTZs20adPH3JychARvv/974f9M2uY7cMPP+Tf/u3fALj++utZuXJlfdsVV1yBiJCb\nm0v37t3Jzc3F4/EwYMAAduzYEfbjnYy3dQz7UtrekhZuF9YZqiIyEXgS8ALPqepDIe0JwIvAUKAY\nuEZVd0Q2qmmND7teyfyDg3j2eB1uPMEKi/EdTtye3OXE7cf7tgZL+6oqEyZM4JVXXmnUZ/369c1e\nyFtPsBbSseV+4eslfzt16sSnn37KsmXLeOaZZ1i0aBEvvPBCi7l8Pl/98YDKyspGj9/SBcaP95zW\nrVvX6ouTH285ZGh8sepjPwOPx9Po5+HxeCK+v/6YgWwhrTIBsAIfS1ocuYuIF3gGuAToD0wVkf4h\n3X4IHFTVbwOPAw9HOqhpnf2Jfcj3xO7ZgyNHjuQf//gHW7fWLktcXl7O5s2b6devH3v37mX16tUA\nlJaWUlNTw5gxY3j55ZcB2Lx5M7t27eLMM8887v0fOHCAYDDIlVdeyQMPPMDatWsBSE1NpbS09Ljf\nl52dzZo1awDqL+MHtfvDn32/tAyHAAAIXklEQVT22fpCWVJS0uT+TvSctm/fXr//OrT4h2vUqFG8\n+mrtvPKXX36Z0aNHt+p+WtLSz+iYF70PMPSrN1rsZ06tcHbLDAe2quo2Va0GXgUmh/SZDCys+/oN\nYLy0dohiIur08gLOC651OsZxZWRksGDBgvqpfyNHjmTTpk3Ex8fz2muvcdtttzFo0CAmTJhAZWUl\nt9xyC4FAgNzcXK655hoWLFjQaIQaas+ePYwdO5bBgwczffp0HnzwQaD24OPNN99cf0A11B133MHv\nfvc7Ro0aVX8dV4CbbrqJXr16MXDgQAYNGlS/5s2MGTO45JJLGDdu3HGfU2JiIvPmzeOyyy5j9OjR\n9O7du1U/s6eeeor58+czcOBAXnrpJZ58MjqrMV577bU8+uijDBkyxA6otkEtLvkrIlcBE1X1prrt\n64ERqnprgz7/qutTWLf9RV2fA83dJ7R+yd9Fq3eT9O4d5AYan5izz9OdX3W4F4A7Kx4nJ9D4AhU7\nPL35dYfaGRT3lj9IVrCwUXuBty+PJd0OwK+P3kc3LWrU/qk3l/9MuhmA3xy9mzRtPKL52DeM5xJv\nBOCZsp8RT1Wj9r/Hjea/EqYiGmTe0duaPK9lcRfyRsJ3SdJynjo6q0n72/GX8Zf4S+kYPMij5f/R\n9OcSP4X/jh/PacF9PFD+QP3tnYMlVBNH1/t2Am13yV8Te469lipmZ1Al8RzydOLXSbPY4c1mjH8l\n11c1/WTyqw6/Yp/nNCZU/w/fq246W2hWh//LIU8nLq9eyuTqprv8bkt+jEpJ4qqqN7nY/16T9hnJ\nT6Pi4ftVr3CBf2WjtmoS+ElK7TkfN1XOZ0RN4/pzRFL5eXLtHufbKp5lYOCzRu37JYNfJt8HwB0V\nT3JmYHOj9l2eLB7ocDcAvyx/hOzgzkbtW7zf5pGk2qmt/z4+hysGnd4kfzgiueRvcyPw0L8I4fRB\nRGYAMwB69eoVxkM31bFDHBVpmZRUNZ5TXOHrSk5GCgD+4kxKqhvPVa6KO52crrXtlUVZlNTENWr3\nx2eS06W2vXx/L0oCKY3agwk9yelce9uRr7KpCZY1avcknk5Op9r2Q5qNT6sbtfs69CAnPQXRICX7\n+jR5XgnJ3clJSyE+6KXkq6btHVK6kZOaQnIgQMn+pu0pqd3ISUkhvSadkqKv20voQ3WPc+ja5DuM\niYxP+/yI+AO1hfD0Lp2Ii08hrTyDksNNX6eZGemk+lJILcugpLRpe3a3NI56U+hQ2o2SsqbtZ3RL\nxe9JJOFId0qONm3P6Z6Cigff4R6UlDdur5F4crrXvkc9B0+npLJxe7knpb49UNKTkqrG7/Fyb2dy\nujWsMY2vxVsV16O+xlQdyKTE33jHSMMak57UuP5EQzgj93OB+1T14rrtXwCo6oMN+iyr6/OhiPiA\nfUCGnuDO7WIdzrCRu4kUey05I9yRezj73FcDOSLSR0TigWuBJSF9lgDT6r6+CvjfExV2Y4wx0dXi\nbhlVrRGRW4Fl1E6FfEFVPxeROdROpl8CPA+8JCJbgRJq/wCYGBXOdD5jTsTGbrEvrHnuqroUWBpy\n270Nvq4Ero5sNBMNiYmJFBcX06VLFyvwplVUleLiYhITE52OYk7ALrPXzmRmZlJYWEhRUVHLnY05\njsTERDIzM52OYU7Ains7ExcXR58+TWcZGGPcxdaWMcYYF7LibowxLmTF3RhjXKjFk5ii9sAiRcDO\nFjuevK7AcZc9aCPa+nOw/M5r68+hreeH6D2H3qqa0VInx4p7tIhIfjhnb8Wytv4cLL/z2vpzaOv5\nwfnnYLtljDHGhay4G2OMC7mxuM9zOkAEtPXnYPmd19afQ1vPDw4/B9ftczfGGOPOkbsxxrR7rizu\nIvKAiKwXkXUi8q6ItO6SJw4RkUdFZFPdc3hTRDo6nelkicjVIvK5iARFpM3MehCRiSJSICJbReRu\np/OcLBF5QUT2110drc0RkSwReV9ENta9fm53OtPJEJFEEVklIp/W5b/fsSxu3C0jImmqeqTu638H\n+qvqzQ7HCpuIXETtmvg1IvIwgKre5XCskyIiZwFB4PfAHaoa81dmqbsY/GZgAlBI7bUMpqrqBkeD\nnQQRGQOUAS+q6tlO5zlZItID6KGqa0UkFVgDfKet/A7qrh2drKplIhIHrARuV9WPTnUWV47cjxX2\nOsk0c8m/WKaq76pqTd3mR0CbW35PVTeqaoHTOU5SOBeDj2mqupzaayq0Sar6paqurfu6FNgI9HQ2\nVfi01rHr88XV/XOk/riyuAOIyK9FZDdwHXBvS/1j2A+Avzkdop3oCexusF1IGyosbiMi2cAQ4GNn\nk5wcEfGKyDpgP/DfqupI/jZb3EXkPRH5VzP/JgOo6i9VNQt4GbjV2bRNtZS/rs8vgRpqn0PMCec5\ntDFhXejdRJ+IpACLgf8T8kk85qlqQFUHU/uJe7iIOLJ7rM2u566qF4bZ9Y/AX4HZUYxz0lrKLyLT\ngMuB8bF6PdqT+B20FYVAVoPtTGCvQ1narbp91YuBl1X1T07naS1VPSQiHwATgVN+gLvNjtxPRERy\nGmxOAjY5laU1RGQicBcwSVXLnc7TjoRzMXgTRXUHJJ8HNqrqb53Oc7JEJOPY7DYRSQIuxKH649bZ\nMouBM6mdrbETuFlV9zibKnx1FxpPAIrrbvqoLc32ARCR7wJPAxnAIWCdql7sbKqWicilwBN8fTH4\nXzsc6aSIyCvAWGpXJPwKmK2qzzsa6iSIyGhgBfAZte9fgHvqruMc80RkILCQ2tePB1ikqnMcyeLG\n4m6MMe2dK3fLGGNMe2fF3RhjXMiKuzHGuJAVd2OMcSEr7sYY40JW3I0xxoWsuBtjjAtZcTfGGBf6\n/4FQNkaHf+z0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1f42e86d6a0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#using a DFT a discrete fourier transform; I think this is the same thing as a fourier series #\n",
    "# discrete sample of the function of interest\n",
    "\n",
    "#create a heaviside function.\n",
    "heaviside_discrete = list();\n",
    "for i in x:\n",
    "    if(i > L/2 or i < -L/2):\n",
    "        heaviside_discrete.append(0);\n",
    "    else:\n",
    "        heaviside_discrete.append(1);\n",
    "\n",
    "dft_heaviside= np.fft.fft(heaviside_discrete); #these should be all the fourier components of the step\n",
    "#to check, let's reconstruct the fourier series\n",
    "max_N = int(len(x)/2);\n",
    "\n",
    "fourier_rec = 0.5/(2*L)+cmath.sqrt(-1)*0.0; #0th order is just the average of the step;\n",
    "for i in range(0, max_N):\n",
    "    index = i+max_N;\n",
    "    fourier_rec += (1/(2*L))*dft_heaviside[index]*np.exp(-cmath.sqrt(-1)*2*np.pi*i*x/(2*L))/max_N;\n",
    "plt.figure()\n",
    "plt.plot(np.real(np.fft.fftshift(dft_heaviside))[:]);\n",
    "plt.title('Fourier components from numpy fftshift')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.plot(x, heaviside_discrete);\n",
    "plt.plot(x, np.abs(np.fft.ifft(dft_heaviside)), '--'); #it's close, but scaling seems off\n",
    "plt.legend(('original', 'reconstructed from fft'))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dft(data):\n",
    "    '''\n",
    "    data is a 1D discrete set of numbers (can be complex)\n",
    "    essential idea is a dft maps one nx1 vector to another nx1 vector, both in the complex plane\n",
    "    basic algorithm as is is O(n^2), which is sub-optimal\n",
    "    :param data:\n",
    "    :return:\n",
    "\n",
    "    '''\n",
    "    i = cmath.sqrt(-1);\n",
    "    num_ord = len(data);\n",
    "    xk = [];\n",
    "    for n in range(num_ord): #iterate through all values of data\n",
    "        sum = 0;\n",
    "        for k in range(num_ord): #for each data point, iterate through all fourier orders\n",
    "            sum+=data[n]*np.exp(-2*np.pi*i*n*k/num_ord)\n",
    "        xk.append(sum);\n",
    "    return xk;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
